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Abstract Numerical MHD simulations of 3D reconnection events in the solar 

corona have improved enormously over the last few years, not only in resolution, 

but also in their complexity, enabling more and more realistic modeling. Various 

f^ ' ways to obtain the initial magnetic field, different forms of solar atmospheric 

O^l , models as well as diverse driving speeds and patterns have been employed. This 

^ ' study considers differences between simulations with stratified and non-stratified 

O ! solar atmospheres, addresses the influence of the driving speed on the plasma flow 

and energetics, and provides quantitative formulas for mapping electric fields and 

dissipation levels obtained in numerical simulations to the corresponding solar 

quantities. The simulations start out from a potential magnetic field containing a 

null-point, obtained from a Solar and Heliospheric Observatory (SOHO) magne- 

Ph ' togram extrapolation approximately 8 hours before a C-class flare was observed. 

C/^ I The magnetic field is stressed with a boundary motion pattern similar to — al- 

(~| ' though simpler than — horizontal motions observed by SOHO during the period 

CIh! preceding the flare. The general behavior is nearly independent of the driving 

speed, and is also very similar in stratified and non-stratified models, provided 

Vh ! only that the boundary motions are slow enough. The boundary motions cause a 

C/3 ' build-up of current sheets, mainly in the fan-plane of the magnetic null-point, but 

^ I do not result in a flare-like energy release. The additional free energy required for 

the flare could have been partly present in non-potential form in the initial state, 

fSJ I with subsequent additions from magnetic flux emergence or from components of 

^ ' the boundary motion that were not represented by the idealized driving pattern. 

00 ■ 
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f^ ■ 1. Introduction 

o 



There have been different attempts to initialize the magnetic field of the pho- 
tosphere and corona for numerical simulations; amongst others by elimination 
of the complex observed small scal e structure by the use of several photo- 



^ I spheric magnetic mono pole sources (IPriest. Buneev. and Titovl . 119971 ). by flux 

^ . emergence ex periments ( Archontis et g/.l . l2004r ) . as well as by extrapolation {e.g. 

'•^ ' A/Taaann ^f nl 9000]! nf anlftr r^KcorTraf r^r-\r TnQcrnAtno'Tiiima ^ n frr»m f^OPTO T^tio 



Masson et al . 2009J) of solar observatory magnetograms, e.g. from SOHO. The 



latter type has typically been used together with potential extrapolations, for 
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simplicity reasons as well as due to the limited availability of vector magne- 
tograms. As potential magnetic fields contain no free magnetic energy, these 
cannot directly be used for explaining how flare events take place and where 
the released energy arises from. Therefore, to use a potential magnetic field as 
the basis for an investigation of a flare event, the field must be stressed into a 
state where it contains sufficient free magnetic energy to account for the energy 
release event. There are different ways by which this may be accomplished. A 
simple approach is to impose boundary mo t ions that resemble the ones de rived 
from observations ( Bingert and Peteil I2OIII : iGudiksen and Nordlundl . 12002 ) . An 



alternative, more challenging approach is to stress the system by allowing ad- 
ditional magnetic flux to enter through photospheric magnetic flux emergence 
(JFan and Gibsonl . 120031 ). In the solar context both of these processes take place 



simultaneously, while experiments typically concentrate on a single type of stress- 
ing, in order to investigate in detail its influence on the dynamical evolution of 
the magnetic field. 

Th e present investigation is an extension of the work done by iMasson et al. 



( 20091 ) . They studied the evolution preceding a specific flare event observed with 



SOHO, starting by taking a magnetogram from about 8 hours before the flare and 
deriving a potential magnetic field. Due to the presence of a 'parasitic' magnetic 
polarity, the resulting magnetic field contains a magnetic null-point. From the 
motions of observed magnetic fragments a schematic photospheric velocity flow 
was constructed, and was used to stress this initially potential magnetic field. The 
imposed stress distorts the magnetic field, causing electric currents to build up 
in the vicinity of the magnetic null-point. The magnetic dissipation associated 
with the electric current allows a continuous reconnection to take place. The 
boundary driving together with the reconnection causes the null-point to move. 
The locations of the magnetic dissipation agree qualitatively with the locations 
of flare emission in various wavelength bands, which may be seen as evidence 
supporting a close association between reconnection at the magnetic null-point 
and th e observed C-c l ass fla re. 

The iMasson et al\ ( 20091 ) paper raises many interesting questions, some of 



which we attempt to answer in the present paper. We therefore let the same 
observations provide the basis for deriving a potential initial magnetic field, and 
employ the same imposed boundary stressing of the magnetic field, using the 
setup to investigate the impact of varying the amplitude of the driving speed, as 
well as the impact of allowing the experiment to take place in a gravitationally 
stratified setting. 

The paper is organized as follows: In Section [5] we list the equations we solve, 
and briefly describe the numerical methods used to solve them. In Section [3] we 
give an overview of the different numerical experiments, in Section U] we present 
and discuss the results, and finally in Section [5] we summarize the main results 
and conclusions. 



2. Methods 

The simulations have be en performed using the ful l y 3D resistive and cornpress - 
ible Stagger MHD code (JNordlund and GalsgaardL Il997at iKritsuk e^~aZI . I2OIII) 
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The following form of the resistive MHD equations 
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solved in the code: 
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where p is the mass density, u the bulk velocity, p the pressure, j the current 
density, B the magnetic field, g the acceleration of gravity, e the thermal energy 
per unit volume and 77 the resistivity. Sij is the shear tensor, r^j the viscous stress 
tensor and fg is a weak diffusive flux of thermal energy needed for numerical 
stability. The term Qy represents viscous dissipation, turning kinetic energy 
into heat, while Qj is the Joule dissipation, responsible for converting magnetic 
energy into heat. 

The solution to the MHD equatio ns is advanced in time using an explicit 3rd 
order predictor-corrector procedure ( Hvmanl ■ ll979^ ■ 



The version of the Stagger MHD code used here assumes an ideal gas law and 
includes no radiative cooling and heat conduction. The variables are located on 
different staggered grids, which allows for conservation of various quantities to 
machine precision. The staggering of variables has been chosen so V • B is among 
the quantities conserved to machine precision. Interpolation of variables between 
different staggered grids is handled by using 5th order interpolation. In a similar 
way spatial derivatives are computed using expressions accurate to 6th order. 
To minimize the influence of numerical diffusion dedicated operators are used 
for calculating both viscosity and resistivity. The viscosity is given by 

v^ M{viCf + V2\n\ +v^M\ - V • u| + ) , (14) 
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where Ad is the mesh size and vi = 0.005 — 0.02, V2 = 0.005 — 0.02, and 
vj, = 0.2 — 0.4 are dimensionless coefficients that provide a suitable amount of 
dissipation of fast mode waves {vi), advective motions {vt), and shocks (fc's). 
The expression | — V • u|+ denotes the positive part of the rate of compression 
—V • u. Cf is the fast mode speed defined by c/ = ^JiB"^ + ^p)/ p- 

The resulting grid Reynolds numbers, A(ic//t^, are on the order of 50 - 200 
in regions with smooth variations, while in the neighborhood of shocks they are 
of the order of a few. The corresponding expression for the resistivity is 

■q = /\d{viCf + V2\vi\ +i/3Ad| - V-u^| + ), (15) 

where uj^ is the component of the velocity perpendicular to B, and where the ex- 
pression scaled by u^^ prevents electric current sheets from becoming numerically 
unresolved. The resulting magnetic grid Reynolds numbers are of the order of a 
few in current sheets, as required to keep such structures marginally resolved. 

The overall scaling with Ad ensures that advection patterns, waves, shocks 
and current sheets remain resolved by a few grids, independent of the mesh size. 

The advantage of these three-part expressions for the viscosity and the re- 
sistivity, compared to having constant viscosity and resistivity is that constant 
values would have to be chosen on the order of the largest of these three term, 
in order to handle shocks and current sheets. In the rest of the volume the 
viscosity and resistiv ity would then be or ders of magnitude larger than needed. 



As demonstrated in (jKritsuk et al\ . 120111) the results are quite similar to state 



of the art codes that use local Riemann solvers. Such codes also have dissipative 
behavior on the scale of individual cells - no numerical code is 'ideal' in the 
sense that it presents solutions corresponding to zero resistivity. 

In this article we refer to the — u x B term in the induction equation as the 
advective electric field, while its counterpart ryj is referred to as the diffusive 
electric field. 



3. Simulations 



The experimental setup is inspired by the work bv iMasson et al\ ( 20091) . Our 



study sets out from a Fast Fourier Transform potential extr apolation applied to a 



level 1.8 SOaO / Michelson Doppler Imager magnetogram (IScherrer et aZ.1 . 119951 ) 
from November 16, 2002 at 06:27 UT, 8 hour prior to a C-class flare occurrence in 
the AR10191 active region. The extrapolation le ads to a 3D magnetic null-point 



topol ogy with a clear fan and spine structure (JGreenl . 1 19891 : IPriest and Titov 



19961 ). In order to allow periodic boundary conditions in the potential field 
extrapolation we applied a windowing function to the SOHO cutout of the active 
region AR10191, which decreases the field close to the boundary towards zero. 
This cut out from the c o mplet e solar disk SOHO data differs slightly from the one 
used bv lMasson et al\ ( 20091) . As a result, in our case the null-point is initially 



located at a height of about 4 Mm above the magnetogram, while in their case 
the null-point is located at a height of only 1.5 Mm above the magnetogram. The 
change in null position places our null-point in the corona proper and allows us 
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Table 1. Simulation runs. In the stratified runs min and max values are given in brackets. 
'Boundary' refers to the plasma flow boundary condition at the lower boundary. Runs with 
a stratified atmosphere are denoted with an 'S', while 'O' stands for open boundary. 



Run 


boundary 


max. driving [kms ^] 




density [cm 


-'] 


temperature [K] 


1 


closed 


3.33 




6.8x10^2 




5.0x105 


lO 


open 


3.33 




6.8x10^2 




S.OxIQS 


2 


closed 


6.67 




6.8x10^2 




5.0x10^ 


3 


closed 


10 




6.8x10^2 




5.0x10^ 


4 


closed 


20 




6.8x10^2 




5.0x10^ 


IS 


closed 


3.33 


[4. 


5x10^0, 9.1X 


10^5] 


[8000,1x10*5] 


3S 


closed 


10 


[4. 


5x10^0, 9.1X 


1015] 


[8000,1x106] 


4S 


closed 


20 


[4. 


5x10^0, 9.1X 


10i5] 


[8000,1 xlO^] 



to perform affordable simulations with stratified atmospheres, while at the same 
time it does not influence the nature of the current sheet formation, and still 
successfully describes a typical solar-like magnetic field geometry. 

Numerically, a magnetic field derived from a potential extrapolation is not 
necessarily divergence free, and we therefore initially apply a divergence cleaning 
procedure, which removes the divergence (as measured by our specific numer- 
ical stretched mesh derivative operators) by iteratively applying a correction 
obtained from solving the Poisson equation 

A(5$ = -V-B. (16) 

The structure of the resulting initial magnetic field is illustrated in Figure [T^ 
and[T|D, showing a strong overall magnetic field and a weaker fan-spine structure 
close to the photospheric boundary, which is separately illustrated in Figure 
m Note, that the field lines of the fan-spine topology were selected specifically 
to show the topology, and that the density of the field lines is therefore not 
representative of the magnetic flux density. 

We performed two types of simulation; one type in which we imposed a 1-D 
gravitationally stratified atmosphere profile, and a second type with constant 
density and temperature. The density and temperature profiles are shown in 
Figure |3l and a summary of the simulation runs may be found in Table [1] 

The magnetic fields are in all cases anchored at the vertical boundaries, which 
due to boundary conditions also prevent plasma from flowing in and out. The 
exception to this is run 10, in which instead constant pressure is assumed at 
the lower boundary, and plasma flows through the boundary are allowed. Only 
minor differences were found between the open and closed boundary cases. 



An imposed horizontal velocity field is introduced at the lower boundary of 
the computational box. This schematic velocity field is based on the motions of 
magnetic fragments observed by SOHO in the active region. These fragments 
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(a) 




(b) 

Figure 1. Magnetic field resulting from the FFT extrapolation of the SOHO magnetogram 
taken on November 16, 2002. The region shown here is the entire computational box, having 
an extent of 60 X 175 X 100 Mm. The slice represents the vertical component of the magnetic 
field. Black is the negative polarity, white is positive. 
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outer spine 



Figure 2. Zoom in of Figurc[T] (upper drawing) to the fan-spine topology, excluding the large 
scale field. The volume below the fan- plane is referred to as 'the dome'. 





10 20 30 40 50 60 
height [Mm] 

Figure 3. Mass densities are presented in black, temperatures in gray. The solid line shows 
the f D mass density and temperature profiles as a function of height x [Mm] of all stratified 
runs. The dashed line illustrates the constant mass density and temperature of run f , run f O, 
run 2 run 3 and run 4. 



outline the fast relative motions obs erved prior to the flare in the active region 
as discussed bv lMasson et al. 



(j2009f) . We used and implemented the description 
of the velocity pattern provided in this reference. In order not to produce initial 
transients the velocity field at the bottom boundary is slowly ramped up, over a 
period of about 100 s, by using a hyperbolic tangent function, and is afterwards 
kept constant. A variety of driv ing speeds ha. v e been employed, ranging from 
values similar to those used bv iMasson et al\ ( 2003 ) to values about 6 times 
lower. We compare and discuss their influence in Section |4l 

In general, it is important to keep the driver velocity well below the Alfven 
velocity of the magnetic concentrations, because the magnetic structure and the 
plasma need to have enough time to adapt to the changing positions of the 
magnetic field lines at the boundary. Ideally, to allow gas pressure to equalize 
along magnetic fields, in response to compressions and expansions imposed by 
the boundary motions, the driver velocity should also be small compared to the 
sound speed in the coronal part of the model. This condition is generally fulfilled 
in all experiments, since the coronal sound speed is on the order of lOOkms"^, 
while our driving speeds are considerably smaller than that. 

With these conditions we ensure an almost force free state at all times in 
the simulation, which implies that the electric current is well aligned with the 
magnetic field. Nevertheless, the line-tied motions of magnetic field lines imposed 
by the lower boundary motions causes the creation of a current sheet in which 
magnetic reconnection takes place. 
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Figure 4. Cell sizes in km for all three axis plotted against the grid indices. X is the height. 



The maximal velocity of around 20kms^^ that we applied exceeds the actual 
velocities measured in the active region by a factor of about 40, while it is at 
the sam e time clearly sub-A lfvenic. This speed up, which is similar to the one 
used bv lMasson et al\ ( 20091 ) . has the desirable effect that we can cover a larger 
solar time interval; a simulated time interval of 12 minutes then corresponds to 
8 hours of real solar time. In our slowest cases (IS, 1, and 10), the simulated 
time is more than an hour, and the driving speed (3.33 kms~^) is approaching 
realistic solar values. 

The simulated region has a size of 62 x 175 x 100 Mm, where our x axis points 
in the direction normal to the solar surface. The computational box is covered 
by a stretched grid of dimensions 320 x 896 x 512, with a minimum cell size 
of sa 80 km maintained in a relatively large region around the null-point. The 
grid size is smaller than 85 km over a 8 x5G x 30 Mm region, which includes the 
entire fan-plane and its intersection with the lower boundary. The distributions 
of cell sizes over grid indices are illustrated in Figure 21 In the initial setup, the 
null-point is located at height index x = 50. 



4. Results and Discussions 

As mentioned above, the field extrapolation based on the SOHO magnctogram 
leads to a fan-spine topology of the magnetic field, illustrated in Figured This 
structure is surrounded by a stronger magnetic field, which extends to much 
larger heights into the corona. We concentrate in the present study on the small 
fan-spine structure, which forms as a consequence of a generally positive polar- 
ity in the active region AR10191 hosting a small ('parasitic') negative polarity 
region. The overlying magnetic field lines, including the ones forming the fan- 
plane, are anchored in the photosphere and build together with the spine a 
rather stable magnetic field structure, keeping the plasma from expanding into 
the upper corona. 

We simulate the motion of magnetic field lines located between the large 
scale negative and positive polarities — hence outside the fan-spine structure — 
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which on November 16, 2002 moved a large amount of magnetic flux towards 
the east sidcJ (left hand side in Figure [T]) of the dome, which we define as the 
volume confined by the fan-plane. This translational motion at the photospheric 
boundary is represented in our experiment by a boundary motion ('driver'), 
which is applied at the lower boundary of our computational box. The boundary 
motions lead to an eastward directed motion of the magnetic field lines outside 
the dome and to magnetic plasma being pushed against the west periphery of the 
fan-spine structure. Especially in the stratified case, a part of this fiow extends 
upward along the magnetic field lines toward the neighborhood of the null-point 
and the outer spine. 

The displacement of field lines, particularly outside of the fan-plane, intro- 
duces a misalignment between the inner and outer spine (see also Figure [5|) . 
Figure [S] shows the field line shear at a nominal 'boundary displacement' (simu- 
lation time times the amplitude of the average applied boundary velocity), D = 
2.15 Mm after the start of run 3, which has a driving speed of about lOkms^^. 
Choosing the displacement instead of the simulation time has the advantage that 
at a given displacement, all runs have experienced about the same energy input 
from the work introduced by the boundary driving and are hence comparable. 
The angle 0, designating the difference of the direction between the inner mag- 
netic field lines (in dark blue) and the outer magnetic field lines (in orange) with 
respect to the fan-plane, is still quite small. The quasi-transparent slices show 
the bulk speed; which is high just outside the fan-plane, where plasma is pushed 
up by the driver. Magnetic field lines closely approaching the null-point run just 
below this high bulk fiow layer. 

The applied photospheric driving motion indirectly moves the fan-plane foot 
points at the west side of the fan-spine structure, causing a slight shear between 
the inner and the outer field lines of the fan-plane to arise due to a different 
stress level of these two fiux systems. The magnetic flux system inside the dome 
experiences a compression of about 5 times the surrounding gas pressure when 
the magnetic flux system west to the dome has moved towards it. This leads to a 
large stress close to the null-point, on the east side of the outer spine, where the 
magnetic field as a consequence reconnects with the surrounding fleld in order to 
reach a lower energy state. A thin current sheet forms in the fan-plane, with the 
largest electric current densities occurring closest to the driver, where the shear 
of the field lines is largest. The magnitude of the electric current is lower in the 
neighborhood of the null-point. At the null-point itself the effect is to disrupt the 
structure of the null in such a way t hat the two spine axes move apart, as seen also 
in other single null investiga tions (jPontin. Bhattachariee. and GalseaardL 12007 : 



Galsgaard and Pontinl 120111) . Figure [S] shows the streamlines and direction of 
the highest electric current in run 3 at D = 2.68 Mm. We find that the electric 
current is mostly anti-parallel to the magnetic field lines in the fan-plane. 

A partial outcome of the reconnection is seen in the motion of the inner 
spine, which is not being moved directly by the applied photospheric driver, but 
nevertheless moves at the photospheric level a significant absolute distance of 



^We use as reference system the solar coordinate system. 
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|u| [krrrsi] 




Figure 5. Magnetic field lines connecting to the inside (dark blue) and to the outside (orange) 
of the fan-plane. The bottom slice shows the bulk speed at a height of about 1.4 Mm. The 
snapshot is taken at a displacement D = 2.15 Mm in the experiment run 3. The box size is 
10 X 16 X 22 Mm. 



[mA'm'^ 




Figure 6. The electric current (red-yellow) flowing antiparallel to the magnetic fleld lines 
(black-gray). The color gradient along the streamlines indicates the flow direction (going from 
red to yellow) and the magnetic field orientation (going from black to gray) . The bottom slices 
show the electric current density of run 3. D = 2.68 Mm and the box size is 10 X 16 X 22 Mm. 
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Figure 7. Relative motion of inner and outer spine for run 3. The top panel shows the cutout 
region for D = 0.03 Mm and the lower panel shows approximately the same selected region for 
D = 0.03, 0.55, 1.08, 1.62, 2.15 Mm. The yellow volume indicates the highest electric current 
density, which is located in the fan-plane. 



about 4.5 Mm in the simulation. The motion is nearly linear in space and time. A 
second signature is the character of the displacement of the outer spine relative 
to the position of the inner spine. 

The initial null-point area, connecting the inner and outer spine, stretches 
with increasing displacement into a 'weak field region', where the magnetic field 
strength is very low. This initiates an electric current that passes through the 
fan-plane, cau s ing th e ratio of the smallest to largest fan-eigenvalues to decrease 



( Parnell et al\ . \l99S \. In this case the null almost adopts a 2D structure. This 



is illustrated in Figure [T] We discuss the relation of this inner and outer spine 
distance to the electric field at the end of section 14.11 

4.1. Time evolution of the diffusive electric field 



The diffusive part of the electric field (yyj) included in the induction equation is 
responsible for both changing the magnetic field topology and for transforming 
magnetic energy into Joule dissipation in the MHD picture. In the Sun the 
diffusive electric field co mponent parallel to the magnetic field is responsible 
for particle acceleration ( Arzner and Vlahosl . l2006r ). It is therefore of particular 
interest to see how this field evolves with the boundary displacement, and to find 
out where it concentrates and how large values it reaches. Figure |8] illustrates the 
electric current accumulation in the fan-plane and along the spine axes of the 
magnetic null. This behavior is representative for all stratified and non-stratified 
runs. In the code the resistivity is not a simple constant, as specified by Equation 
(fT4)) . allowing diffusion to be locally increased where dissipation is needed to keep 
structures from becoming unresolved, while at the same time allowing a minimal 
amount of diffusion in regions where the magnetic field is smooth. Images of the 
diffusive electric field are therefore not exact replicas of images of the electric 
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Figure 8. Current density in y-z slice for constant density run 4 (upper row) and stratified 
atmosphere run 4S (lower row) at a height of x = 3.22 Mm and different displacements D. 



current, but since the resistivity is generally near its largest value in current 
sheets there is a close correspondence. 

The diffusive electric field is found to be concentrated in the fan-plane and 
to have its local peak in the region where the fan-spine intersection is distorted. 
However, large values occur over a significant fraction of the fan-plane, as is the 
case for the electric current density. 

It is only the parallel diffusive electric field that gives rise to both magnetic 
reconnection and particl e acceleration and its magn i tude indicates how violent 
these processes can be ( Schindler. Hesse, and Birnl 119881) . In our simulations 
the advective electric field (— u x B), associated with the bulk plasma motion 
in the fan-plane and along the spine axis is much stronger than its diffusive 
counterpart (Figure [5a|) . but since it is perpendicular to the magnetic field 
this component causes no magnetic dissipation, and cannot be associated with 
particle acceleration. 

When determining the values of the diffusive electric field one finds that the 
peak values in the vicinity of the null-point are increasing with growing dis- 
placement (see Figure l9b| . going from initially zero to on the order of 50 Vm~^, 
while in the fan-plane the diffusive electric field reaches more than 90Vm~^ 
at D = 2.95 Mm. However, as shown below (c/. Equation [^ and Equation [531) , 
to estimate the analogous solar electric field, the simulation value should be 
reduced with a factor equal to the power of 0.3 of the factor (about 20 for run 
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(b) 



Figure 9. Maximum advcctive electric field (a) and average diffusive electric field calculated 
over those grid points at which the electric field amplitude is within 10 % of the maximum 
field in a slice of 2 grid cells thickness (b) of the cutout region at the respective height of the 
null-point (*) and in the total cutout region 10 X 16 X 22 Mm around the null-point, including 
the fan-plane (+) in run 3. 



3) by which the boundary driving is exaggerated; here we obtain E w 90/20°'^, 
or a bout 36Vm~^. Consid ering our more benign conditions, this is consistent 
with lPudovkin et all (1998|), who find typical electric fields to be of the order of 
100 - 300 Vm"^ under flaring conditions. Electrons accelerated along the entire 
current sheet, with an extent of about 15 Mm (see Figure |6]), could nevertheless 
gain energies of up to about 300 MeV in our case. Of interest is also the shape 
of the average diffusive electric field increase in the vicinity of the null-point, 
plotted in Figure |9b] as a (*)-line. Its progression shows an almost identical 
behavior as a plot of the increasing distance between the inner and outer spine 
(see Figure [7]) plotted against the displacement (plot not presented here). 



4.2. Comparison of stratified and non-stratified simulations 

Since the null-point is an essential node of the magnetic skeleton, we use it as 
a reference point for our investigation of the influence of the density profile on 
the temporal changes of the magnetic field. Here we compare simulation run 4S, 
which has a stratified atmosphere, with simulation run 4, which has a constant 
chromosphere- like density and temperature atmosphere (see. Table [!}. Figure [TU] 
shows the position of the null-point, connecting the fan-plane and spine magnetic 
field lines, in all three directions versus the boundary displacement D. In the 
stratified case the very dense plasma (9 x 10^^ cm~'^) at the bottom of the box 
gives rise to a low Alfven speed (v^ = B/y^/ToPm); meaning that the higher 
density impedes the propagation of the boundary disturbance into the box. But 
once the disturbances reach beyond the transition region the ambient density 
has fallen drastically, and the dense plasma from the lower atmosphere starts 
spreading out faster and pushes the null-point up. The velocity gradient induced 
by the large density change causes strong dissipation and in some regions a 
disruption into multiple smaller current sheets, which can also be seen in the 
electric current density comparison in Figure [51 showing \j\ for runs 4 and 4S. 
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Figure 10. Null-point position in Mm in the coordinate system of the experiment. The vertical 
coordinate is x, and the two horizontal coordinates are y and z. The solid line shows the 
constant density run 4, while the dashed line represents the stratified run 4S. 



In the non-stratified runs the amount of plasma that moves upwards along 
the magnetic field lines into the corona and contributes to pushing the null-point 
upwards is much reduced, as seen in Figure [TU] 

For all runs we find a significant null-point motion (about 1-2 Mm in x, 6 Mm 
in y and 1 Mm in z) due to the applied driving motion on the bottom boundary, 
which is comparable to the relative motion of the inner spine. 

4.3. The influence of boundary conditions 

The driver of the magnetic field evolution is the boundary motion. The magnetic 
field displacements imposed by the boundary motions have a large influence on 
the spatial structure of the magnetic skeleton, consisting of null-points, sepa- 
ratri x surfaces (such as the fan-plane), se parators, sources and flux domains 
fe.q. lParnell. Havnes. and GalsgaardluOOSl and references therein). The skeleton 
itself is a very robust structure, which does not change from a topological point 
of view, but the detailed appearance of it changes with boundary displacement, 
as already shown by the analysis of the null-point motion. 

In Figure fTTbl we compare the density profiles of the closed runs 1,2,3 and 4 
and the open boundary run 10. The figure shows that all non-stratified closed 
fiow boundary runs develop a similar density profile: A certain expansion or 
compression of a region connected to the closed boundary leads to the same 
density profile, regardless of the driving speed. 

The initial sound speed is about 83kms~^ for the non-stratified runs, while 
being much lower — on the order of lOkms"^ — in the lower parts of the stratified 
runs. The sound speed influences among other things the in- and outflows at open 
boundaries. 



Figure [llal illustrates that the closed boundary conditions influence particu- 
larly the lowest density layer, where regions with low and high density build up 
and cannot be emptied nor flUed by plasma out- and inflow. This is not very 
important for the dynamics of the system, but is also not very solar-like. In the 
open case, there is a continuous mass exchange and pressure equalization at the 
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(a) 



(b) 



Figure 11. Density profile min and max over each slice in height for the closed non-stratified 
simulations, normalized to the asymptote: run 1, run 2, run 3 and run 4 and the open non-s- 
tratified run lO. (a) is a zoom in of the first 0.5 Mm of plot (b). The displacement of each run 
is taken to be from the snapshot closest to D = 1.05 Mm. 



boundaries, in which case it is crucial that the sound speed is well above the 
driving speed, so that the system has time to approach pressure balance. 

The low sound speed in the stratified case poses a very tight restriction on 
the driving speed, in order to avoid exaggerating the effects of inertia. On the 
other hand there is a clear advantage of having stratification: it provides a pool 
of mass for the corona to communicate with; the large amount of mass at low 
temperature acts as a buffer, due to the low Alfvcn and sound speed. 



Overall, Figure lllbl confirms that the density contrast is mainly caused by 
volume changes, which arise from the imposed boundary motions. If these volume 
changes happen sufficiently slowly relative to the Alfvcn and sound speed, the 
driving speed loses its importance for the results (but not for the computational 
cost of obtaining them!). 



4.4. Energy dissipation 

As the boundary moves according to the prescribed driving pattern, with mag- 
netic field lines passing through the boundary essentially 'frozen in', because 
of the boundary conditions, the system response may be mainly split into two 
distinct components. The first is the change in potential magnetic field energy 
due to the change of the vertical magnetic field component brought about by the 
boundary motions; this is the smallest amount by which the magnetic field energy 
could change. Secondly, in addition, the 'free magnetic energy' component will 
change as well. This non-potential part of the magnetic field is (by definition) 
associated with a non-zero electric current proportional to V x B. The electric 
current may either be smooth and space-filling, or may be concentrated in electric 
current sheets, corresponding to near-discontinuities of the magnetic field. 

Formally, the rate of change of magnetic energy density cb = B^ /2 is de- 
scribed by Equation (flT)) . which shows that changes of magnetic energy are due 
to the net effect of a (negative) divergence of the Poynting flux Fp, conversion 
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Figure 12. Comparison of the average magnetic energy a) and average magnetic dissipation 
divided by the normalized (vq = 10 km s~^) average boundary driving speed of each run 
b) for constant density runs (sohd hncs) and stratified runs (dashed Unes) in the cutout of 
approximate size 10 x 16 x 22 Mm. 



into bulk kinetic energy by the Lorentz work Wl , and conversion to heat through 
Joule dissipation Qj. 



at 



./• 



(17) 



The Poynting flux is defined as Fp = E x B, and the Lorentz work is Wl = 
u • (j X B). Joule dissipation, Qj ~ E,, • j = Tyj^, primarily takes place in the 
strong current sheets. Electric currents flow mainly along the magnetic field 
in the corona and therefore Qj is a suitable indicator for locations at which a 
significant component of the electric field parallel to the magnetic field may exist. 
Such a parallel electric field can accelerate charged particles along the magnetic 
field lines, resulting e.g. in the brightening of fiare ribbons. 

The evolutions of the magnetic energy and the Joule dissipation normalized to 
the average driving speed are summarized in a plot covering several simulation 
runs in Figure [TH The first thing to notice is that the evolution of the magnetic 
energy, when expressed in terms of the boundary displacement D, is practically 
identical in all of the runs. The reason for this is that most of the boundary 
work goes directly into increasing the potential magnetic energy, while only a 
small amount goes into free magnetic energy. From Figure 112b is seen that 
the dissipation increases with increasing displacement. This indicates that an 
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increasing amount of free energy becomes available through the build up of 
current structures in the null-point fan-plane as the experiment progresses. 

A comparison between different power-law relations indicates, that the nor- 
malization of the magnetic dissipation by the average normalized driving speed 
to the power 0.6 employed in Figure [T2b brings the curves showing the evolution 
of magnetic dissipation for all the different runs closest together into a relatively 
tight set of parallel relations. This illustrates that the rate of magnetic dissipa- 
tion, at any given value of the displacement, is approximately proportional to 
the rate of boundary displacement to the power 0.6. 

The dissipation curves corresponding to the different experiments follow the 
same general trend, although with some differences, in particular between the 
stratified and non-stratified cases. The dissipation is generally higher for the 
stratified cases than for the non-stratified runs during early times and lower 
during late times. The exception is run IS, the run with the lowest driving 
speed, which agrees closely with the non-stratified cases with the slowest driving 
speed. 

The deviations from this common asymptotic behavior are likely consequences 
of the low Alfven speeds in the dense layers of the stratified models, causing the 
dissipation in the stratified runs to be initially high due to their higher densities 
at low heights compared to the constant density cases. Later, when the motions 
introduced by the driver reach greater heights, where the mass density is lower 
than in the constant density runs, the stratified runs generally display a lower 
dissipation. 

We note in this context also that the viscous dissipation in the system is much 
smaller than the Joule dissipation, as is expected in a coronal environment. 

The results summarized in Figure [T^] illustrate that for a quantitatively ac- 
curate estimate of properties related to the magnetic dissipation it is essential 
to drive in a way which is compatible with the ordering of characteristic speeds 
in the Sun; i.e., to keep the boundary speed smaller than the Alfven speed, and 
to scale the quantities down in proportion to the speed-up factor used in the 
driving. 

For the non-stratified runs we compute initial Alfven speeds of 70 - 1400 km s^^| 
at the lower boundary, which is clearly higher than the driving speed of all runs. 
So we expect, as is also shown by Figure [T^ a similar dissipation increase with 
increasing displacement, after some initial differences due to the driving speed 
differences. 

In the stratified atmosphere runs the Alfven speed increases with height. At 
the lower boundary the initial Alfven speed is approximately 5 -^ 200kms~^, 
where the minimum value falls below the driving speed used in runs 3S and 4S. 
Run IS is just at the edge of being driven slower than the minimum Alfven 
speed and indeed gives results which are similar to the non-stratified case run 
1 . Figure [13] shows that the stratified and non-stratified cases are nevertheless 
distinguishable. The volume renderings show the Joule dissipation normalized 
by the driving speed, for run 1 (upper panel) and run IS (lower panel). The 
locations of the dissipation maxima agree nicely, but the dissipation maxima 
differ by a factor of about 1.8, being higher in the stratified run IS. 
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Figure 13. The highest joule dissipation of run 1 in the top row and run IS in the lower 
row is shown as a semitransparcnt volume rendering at about D = 0.88 Mm. The box size 
is 10 X 16 X 22 Mm. The field lines represent magnetic field lines passing closest to the null 
point (purple), in the close surrounding of the null-point (magenta) and the overlying strong 
magnetic field (dark blue). The slice on the lower boundary shows the Bx magnetic field 
component. 



In summary, wc find that the ratio of the driving speed to the Alfven speed 
has a noticeable impact on the dissipation level and, as Figure[T2]illustrates, that 
the stratified simulations tend to display a progressive growth of deviations from 
the common asymptotic relations defined by the non-stratified runs, unless the 
driving speed is small compared to both the local Alfven speed and the sound 
speed. 

However, these are relatively small deviations, compared to the main trend, 
which is a proportionality between the magnetic dissipation and the driving 
speed to the power 0.6. 

4.5. Scaling of magnetic dissipation in the current sheet 

As illustrated by Figure [SJ the magnetic field line orientations on the two sides of 
the fan-plane only differ by a small amount. This means that the electric current 
carried by the current sheet, whose thickness (on the order of a few grid cells in 
the present numerical model) we denote with As, is much less than the maximal 
electric current density 2B/As that would result from a complete reversal of the 
magnetic field orientation across the current sheet. The electric current density 
in the current sheet (CS) is in fact on the order of 



jcs ^ AS As~^ sa sin{(j>)Bcs^s~'^ ~ (f>Bcs^s' 



(18) 



for small 0, where Bqs is a typical strength of the magnetic field just outside 
the current sheet, and </> is an angle characterizing the difference of direction of 
field lines on the two sides of the current sheet. 

A fundamental question is now how the total dissipation in the current sheet, 
and hence the average rate of reconncction in the structure, depends on factors 
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such as the numerical resolution and rate of work done at the boundary. By 
construction the code keeps current sheets just barely resolved. A change in 
numerical resolution is thus directly mapped into a proportional change of the 
current sheet thickness As. This behavior is obtained by making ry essentially 
proportional to the grid size ( cf. Equation llSp . To a first order approximation 
both Bqs and (p are independent of As. By Equation [18] the electric current 
density is inversely proportional to As, and the magnetic dissipation rate per 
unit volume Qj therefore scales as 

Qcs = vfcs - ^B^^s-' = 02b2,^As-i. (19) 

To obtain the total dissipation in the current sheet this needs to be multiplied 
with the volume of the current sheet, 

Fes - AcsAs, (20) 

where we denote by Acs the total area of the current sheet. We thus conclude 
that the total dissipation is 

Qcs^cs - AB^Acs, (21) 

and hence is, to lowest order, independent of As and the resistivity in the current 
sheet. Note that, as a consequence of being small, reconnection in the current 
sheet can proceed without requiring super- Alfvenic outflow velocities from the 
current sheet. 

Estimating now the diffusive electric field in the current sheet we find that it 
scales as 

Ecs = ^/csjcs « AB = 4>Bcs, (22) 

again independent of As, but proportional to the change of magnetic field 
direction across the current sheet and hence proportional to (j). 

Generally the work done by the boundary must go into an increase in magnetic 
energy (potential plus free magnetic energy), or into kinetic energy or ohmic 
dissipation. In the present case the magnetic dissipation is able to nearly keep 
up with the free energy input, and the system essentially goes through a series 
of states not far from potential. As Figure fT2l shows, the total magnetic energy 
depends mainly on the displacement and very little on the driving speed itself, 
while the dissipation is essentially proportional to how fast we drive at the 
boundary to the power of 0.6, thus Qcs oc u°'^. So we conclude with the help of 
Equations ((T9| and ((22)) that the magnitude of the parallel electric field along 
the current sheet scales as 

-Ecs fx </) ex v°-^. (23) 

This electric field - driving speed relation is not an artifact of the chosen nu- 
merical method; the scaling of rj with As is generic to all numerical methods, 
and serves to ensure that higher numerical resolution can be used to reach larger 
magnetic Reynolds and Lundquist numbers. 
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With much higher rates of stressing, or with much higher numerical reso- 
lution, the current sheet may need to fragment and enter a turbulent regime, 
in order to support the required amounts of dissipation and reconnection in 
the face of increasing constraints by the thinness of the current sheets. This is 
a process th at is by now well unders t ood to be able to take over when the 

need arises ( Galseaard a nd Nordhmd. 1996t Nordlund and Galsgaard. 1997b : 

Gudiksen and Nordlund .12005: Bereer aiid Asgari-Targhi . 2009 : Bingert and Peten .| 



20111: iNg. Lin, and Bhattacharied . 120121: iPontinl . 120111 ). 

We thus conclude that, provided that conditions that apply equally to both 
numerical simulations and the Sun are fulfilled, neither the total magnetic dis- 
sipation in a current sheet structure, nor the diffusive part of the electric field 
depend, to lowest order, on the electrical resistivity — or for that matter on the 
precise mechanism that sets the level of the electrical resistivity. 



4.6. Triggering of rapid energy release 



As in lMasson et al\ ( 20091 ) there is no significant and sudden relaxation of the 
system, with an energy release that could correspond to a flare, even though 
the simulations cover enough solar time to get across the observed flaring event. 
Looking at the potential part of the magnetic field at different displacement steps 
in the simulations, we find steadily increasing magnetic potential energies, and 
most of the build-up of magnetic energy seen in the experiment actually goes 
into increasing the potential rather than the non-potential part of the magnetic 
energy. As discussed above, the magnetic dissipation is able to keep up with 
varying levels of boundary work, with a residual amount of free energy scaling, 
if it is proportional to the rate of dissipation, approximately as the driving speed 
raised to 0.6. The stress in the system is demonstrably moderate, even with our 
exaggerated driving speed, and we therefore cannot expect boundary motions of 
the type applied here to be able to explain violent events similar to the observed 
solar flare. 

We consider four possibilities for this difference in behavior between these 
MHD simulations and the Sun: 

i) The limited resolution of the numerical experiment is preventing an instability 
from occurring that would otherwise trigger a flare-like event. 

a) A flare-like event would occur if taking into account kinetic effects {e.g. by 
using particle- in-cell simulations). 

Hi) A flare-like event could take place, with MHD alone, but would require an 
additional Poynting flux through the boundary, in addition to the Poynting: 
flux generated by the simple driver implemented here and in iMasson et al\ 
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iv) The additional free energy available because the system was initially not in a 
potential state could have helped. 

As demonstrated above, numerical resolution does not to a first order deter- 
mine the level of dissipation in current sheets, once they are reasonably well 
resolved. Indeed, by the arguments in the preceding subsection we expect the 
level of stress (as measured for example by the angle 4>) caused by the type of 
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motions we employ in the current investigation to be smaller in the Sun (by 
a factor of e.g. « T^"^ « 2 relative to our runs 1, 10 and IS), and it is thus 
very unlikely that the flare was triggered by accumulation of stress from this 
particular type of boundary motion. 

An MHD-instability could in principle occur at a later point in time than 
to where our runs go; but as demonstrated by our experiments even a driving 
speed that is highly exaggerated relative to the solar value is not able to build up 
sufficient free energy to account for a C-class flare: The maximum total magnetic 
energy in the entire simulated domain of our experiment is on the order of 
2 X 10'^° ergs, of which only a very small fraction is free energy. Estimates o f 



C-flare emission are larger than even our potential energy (jKretzschmaii 120111) , 
and hence the free energy available in our model is not sufficient to power a C- 
class flare. As discussed above, the level of stress is expected to be proportional 
to the driving speed to some small positive power, and is expected to be largely 
independent of the level of resistivity, and hence it is unlikely that more stress 
would build up in the Sun than in these numerical simulations. 

Figure[T31 which is representative of all runs, illustrates that the rate of change 
of magnetic energy increases in the beginning, as a result of the work done by the 
lower boundary. This energy input loses efficiency as the angle between the field 
lines and the driving boundary approaches 90 degrees. Additionally, since the 
driving pattern location is fixed, the applied stress becomes less and less efficient, 
as the flux in this particular region is gradually removed. The dissipation is 
small compared to the rate of change of magnetic energy, implying that the 
dissipated free energy only makes up a tiny fraction compared to the change of 
the potential energy. This is in agreement with Figure [SJ which illustrates that 
the shear of the magnetic field lines in the current sheet is not very large, even 
with our exaggerated driving speed. Hence the reconnection is rapid enough to 
keep the boundary work at a given displacement at an approximately constant 
level, and to keep the system in a near potential state. We are thus far away 
from an instability to occur and it is improbable to achieve one at a later point 
of the experiment with the applied boundary motion. The mostly compressive 
boundary driving pattern probably does not represent reality accurately enough, 
and additional shearing, twisting, or emerging motions may have been present. 

By similar arguments kinetic effects are not expected to have a major influence 
on the macroscopic d issipation b e havio r and reconnection rate. This is consistent 



with the findings by iBirn et all (J2005r ) , which show that the same amounts of 



energy release take place for the same type of reconnection event, independent of 
the dissipation mechanism. It also shows that there are only small differences in 
the time scales of the event going from MHD to PIC simulations. These results 
are entirely consistent with the current sheet scaling arguments presented above, 
since they hold true independent of the nature of the dissipation mechanism. 

We thus are left with the option that the additional free energy mentioned in 
the third and forth alternative above is the most likely explanation for the ob- 
served solar flare. This conclusion is consistent with the fact that all observational 
evidence associates flaring either with new flux emergence, or with magnetic 
configurations that evolve into states that are decidedly MHD- unstable; i.e., ones 
in which lower energy states become accessible and where dissipative effects make 
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Figure 14. The average rate of change of magnetic energy for run 3 is plotted as a solid line, 
and the dissipation is shown as a dashed line. Both are restricted to the cutout of approximate 
size 10 X 16 X 22 Mm used previously. 



it possible to reach those lower energy states and dynamic instabilities therefore 
develop. Some amount of free magnetic energy was clearly present in the system 
at the time where we are forced by circumstances to assume a potential initial 
state. Whether that extra free energy alone would have been enough to cause 
a flare several hours later is an open question. Our driver pattern might have 
been able to add slightly more magnetic energy to a configuration that already 
contained free magnetic energy, but that effect is probably small. We consider 
the addition of Poynting flux, increasing the free magnetic energy of the system 
at a considerably higher rate than can be achieved with our driver pattern as the 
more likely alternative. Additional Poynting flux across the lower boundary can 
take the form of either emerging flux (Poynting flux due to vertical velocity) or 
twisting boundary motions (Poynting flux due to vortical horizontal velocity). 

As evidence for this view, we consider magnetograms starting on November 
12, 2002 in the early afternoon, when a strong negative polarity first appeared 
within the positive polarity. The negative polarity expanded and several erup- 
tions could be seen at EUV wavelengths. During the interval simulated here one 
observes an additional increase in the negative flux inside the dome, and there 
are also other significant rearrangements that our driver pattern is not able to 
represent or explain. 



5. Conclusions 

The main topic of the presented work is a study of the influence of the driving 
speed for stratifled and non-stratified atmosphere models used to simulate 3D 
reconnection events in the solar corona. The major findings are: 

• Sufficiently low driving speeds lead to similar plasma behavior, including 
similar evolution of the magnetic energy density (after compensating for 
the expected effects of the different driving speeds). 
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• The magnetic dissipation and the diffusive electric field — capable of ac- 
celerating charged particles depend only weakly on the grid resolution in 
the numerical experiments. 

• When driving is taking place sufficiently slowly, the rate of magnetic dis- 
sipation increases approximately as the boundary driving speed raised to 
the power 0.6 while the diffusive electric field increases approximately as 
the driving speed raised to the power 0.3. 

• The driving speed has a larger impact on the general plasma behavior in 
the cases with stratified atmospheres, while its influence is minor for non- 
stratifled runs with closed flow boundary conditions, if compared at the 
same boundary motion displacements. 

• The sound speed compared to the driving speed determines the exchange 
of plasma at the boundaries. Voids can be filled in open boundary runs, 
while pressure differences at the boundary are carried along in the closed 
boundary cases. This is not an issue in the stratified runs, due to the large 
reservoirs of relatively cold gas in the vicinity of the lower boundary. On 
the other hand, the lower information speeds there set severe restrictions 
on the driving speed. 

Our investigations suggest that the applied simple driver pattern is unlikely 
to be able to cause a flare-like energy release in the simulations (as well as in 
the Sun). In fact, based on a comparison of the available free energy in our 
model with estimates from observations of C-class flares we anticipate that the 
corresponding solar configuration must have had significantly more free magnetic 
energy added through the boundary, or must otherwise have been in a strongly 
non-potential state already before the studied time interval. 
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